# Required libraries
library(stats)
library(pracma)
# Black-Scholes call price formula
black_scholes_call <- function(S, K, T, r, sigma) {
d1 <- (log(S / K) + (r + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
d2 <- d1 - sigma * sqrt(T)
call_price <- S * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
return(call_price)
}
# Function to find the implied volatility
implied_volatility <- function(S, K, T, r, market_price) {
objective_function <- function(sigma) {
black_scholes_call(S, K, T, r, sigma) - market_price
}
# Use uniroot to find the implied volatility
result <- tryCatch({
uniroot(objective_function, lower = 0.001, upper = 5)$root
}, error = function(e) NA)
return(result)
}
# Given data
S <- 100 # Current index price
r <- 0.03 # Risk-free interest rate
# Market prices matrix
market_prices <- matrix(c(
40.2844, 42.4249, 50.8521, 59.1664,
30.5281, 33.5355, 42.6656, 51.2181,
21.0415, 24.9642, 34.4358, 42.9436,
12.2459, 16.9652, 26.4453, 34.789,
5.2025, 10.1717, 19.4706, 27.8938,
1.3448, 5.4318, 14.4225, 23.3305,
0.2052, 2.7647, 11.2103, 20.7206,
0.0216, 1.4204, 9.1497, 19.1828,
0.0019, 0.7542, 7.741, 18.1858
), nrow = 9, byrow = TRUE)
# Strikes and maturities
strikes <- c(60, 70, 80, 90, 100, 110, 120, 130, 140)
maturities <- c(0.25, 0.5, 1, 1.5)
# Create a matrix to store implied volatilities
implied_vols <- matrix(NA, nrow = length(strikes), ncol = length(maturities),
dimnames = list(paste0("Strike_", strikes), paste0("Maturity_", maturities)))
# Calculate implied volatilities for all combinations
for (i in seq_along(strikes)) {
for (j in seq_along(maturities)) {
implied_vols[i, j] <- implied_volatility(S, strikes[i], maturities[j], r, market_prices[i, j])
}
}
# Print the implied volatilities matrix
print(implied_vols)
## Maturity_0.25 Maturity_0.5 Maturity_1 Maturity_1.5
## Strike_60 NA 0.5865828 0.8175021 0.9478482
## Strike_70 0.2486434 0.5220188 0.7158389 0.8294040
## Strike_80 0.3148182 0.4545265 0.6174672 0.7155847
## Strike_90 0.2801724 0.3894853 0.5277066 0.6147351
## Strike_100 0.2426795 0.3368458 0.4608755 0.5465596
## Strike_110 0.2173596 0.3070502 0.4299974 0.5243019
## Strike_120 0.2056962 0.2975798 0.4272046 0.5335653
## Strike_130 0.2022396 0.2989333 0.4376976 0.5558278
## Strike_140 0.2030054 0.3049178 0.4531255 0.5819529
library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# Prepare the data for plotly
df <- expand.grid(Strike = strikes, Maturity = maturities)
df$Volatility <- as.vector(implied_vols)
plot_ly(data = df, x = ~Maturity, y = ~Strike, z = ~Volatility, type = "mesh3d") %>%
layout(scene = list(
xaxis = list(title = "Maturity"),
yaxis = list(title = "Strike"),
zaxis = list(title = "Implied Volatility")
))